2.2.2 Main effect and post-hoc analysis
# sim_main_posthoc is available in R/funcs.R
p_main_posthoc <- sim_main_posthoc(N_subj = 30, iter = Nsim, n_core=16,
file_cache = here::here("simulation", "p_main_posthoc.rds"),
N_levels = c(3,4,5,6))
2.2.2.1 FWER and CER
df_TypeI_posthoc <- sig_main_posthoc(p_main_posthoc, 0.05) %>%
group_by(N_level, alpha, adjust) %>%
summarize(TypeI_main = mean(sig_main),
TypeI_postonly = mean(sig_anypost),
`TypeI_main&posthoc` = mean(`sig_main&posthoc`),
.groups = "drop")
df_TypeI_posthoc
Claim significant results only when both the main effect and at least one of the post-hoc tests (with different multiple comparison corrections) are significant:
df_TypeI_posthoc %>%
filter(adjust != "uncorrected") %>% # only when using multiple comparison correction
select(-c(TypeI_main, TypeI_postonly)) %>%
pivot_wider(c(N_level, alpha),
names_from=adjust, values_from=`TypeI_main&posthoc`)
Claim significant results only when the main effect and at least one of the post-hoc tests (without multiple comparison corrections) are significant:
df_TypeI_posthoc %>%
filter(adjust == "uncorrected") %>% # only when NOT using multiple comparison correction
select(-c(TypeI_main, TypeI_postonly)) %>%
pivot_wider(c(N_level, alpha),
names_from=adjust, values_from=`TypeI_main&posthoc`)
Claim significant results only when at least one of the post-hoc tests (with different multiple comparison corrections) are significant (regardless of the main effect results):
df_TypeI_posthoc %>%
filter(adjust != "uncorrected") %>% # only when using multiple comparison correction
select(-c(TypeI_main, `TypeI_main&posthoc`)) %>%
pivot_wider(c(N_level, alpha),
names_from=adjust, values_from=TypeI_postonly)
df_posthoc_tmp <- df_TypeI_posthoc %>%
pivot_longer(starts_with("TypeI_"), names_to = "effects", values_to = "TypeI") %>%
mutate(adjust = if_else(effects == "TypeI_main", "uncorrected",
if_else(adjust == "uncorrected", "uncorrected",
str_to_title(as.character(adjust)))),
adjust = fct_reorder(adjust, TypeI, mean, .desc=TRUE),
adjust = fct_relevel(adjust, "uncorrected"),
effects = factor(effects,
levels = c("TypeI_main", "TypeI_postonly", "TypeI_main&posthoc")),
effects = fct_recode(effects, `Main effect\n` ="TypeI_main",
`Pairwise comparisons\n(six tests, FWER)`="TypeI_postonly",
`Main effect & Pairwise comparisons\nwith six tests (CER)`="TypeI_main&posthoc")) %>%
distinct()
fig_simu_b <- df_posthoc_tmp %>%
filter(N_level == 4) %>%
ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) +
geom_col(width=1) +
facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x")+
ggh4x::facetted_pos_scales(
list(
scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
scale_x_continuous(breaks = 3.5, expand = c(0, 0.65)),
scale_x_continuous(breaks = 3.5, expand = c(0, 0.65))
)) +
scale_fill_manual(values = custom_clr,
breaks = levels(df_posthoc_tmp$adjust)) + #,
geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
scale_y_continuous(limits = c(0, 24), expand = c(0, 0)) +
labs(x = "One-way ANOVAs with four levels", y = "False Positive Rate (%)", fill = "") +
theme(legend.position = c(0.74, 0.63),
axis.ticks.length = unit(.15, "cm"),
axis.text.x = element_blank(),
panel.spacing.x = unit(0, "pt"),
strip.placement = "outside",
axis.title.x = element_text(margin = margin(t = 0)))
# ggsave(file.path("figures","fig_simu_TypeI_mainpost.png"), fig_simu_b, width = 8, height = 5.6)
fig_simu_b
for (i in unique(df_posthoc_tmp$N_level)) {
nam <- paste("fig_simu_b_tmp", i, sep = "")
assign(
nam, df_posthoc_tmp %>%
filter(N_level == i) %>%
ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) +
geom_col(width=1) +
facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x")+
ggh4x::facetted_pos_scales(
list(
scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
scale_x_continuous(breaks = 3.5, expand = c(0, 0.65)),
scale_x_continuous(breaks = 3.5, expand = c(0, 0.65))
)) +
scale_fill_manual(values = custom_clr,
breaks = levels(df_posthoc_tmp$adjust)) + #,
geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
scale_y_continuous(limits = c(0, 45), expand = c(0, 0)) +
labs(x = paste0("One-way ANOVAs with ", i, " levels"), y = "False Positive Rate (%)", fill = "") +
theme(legend.position = c(0.75, 0.6),
axis.ticks.length = unit(.15, "cm"),
axis.text.x = element_blank(),
panel.spacing.x = unit(0, "pt"),
strip.placement = "outside",
axis.title.x = element_text(margin = margin(t = 0)))
)
}
fig_simu_b_supp <- ggarrange(fig_simu_b_tmp3, fig_simu_b_tmp4,
fig_simu_b_tmp5, fig_simu_b_tmp6,
nrow = 2, ncol = 2)
# ggsave(file.path("figures","fig_simu_TypeI_mainpost_supp.png"), fig_simu_b_supp, width = 13, height = 8)
fig_simu_b_supp
2.2.2.2 Single tests
# with omnibus F-test and significant effect
df_sig_main_single1 <- sig_main_posthoc(p_main_posthoc, 0.05) %>%
filter(adjust == "uncorrected") %>%
transmute(Method = "Main effect & Pairwise",
N_level, contrast,
sig_cber = `sig_main&posthoc`, # main and post-hoc (uncorrected)
sig_cber_single = sig_cber & sig)
df_TypeI_main_single <- sig_main_posthoc(p_main_posthoc, 0.05) %>%
filter(adjust != "uncorrected") %>%
transmute(Method = str_to_title(adjust),
N_level, contrast,
sig_cber = `sig_anypost`, # post-hoc only (corrected)
sig_cber_single = sig_cber & sig) %>%
# combine two df and make it to long format
bind_rows(df_sig_main_single1) %>%
group_by(N_level, contrast, Method) %>%
summarize(TypeI_cber = mean(sig_cber), # correction and at least one sig
TypeI_cber_single = mean(sig_cber_single),
.groups = "drop") %>%
pivot_wider(c(N_level, Method, TypeI_cber),
names_from = contrast, values_from = TypeI_cber_single) %>%
pivot_longer(!c("N_level", "Method"),
names_to = "contrast", values_to = "TypeI") %>%
filter(!is.na(TypeI))
# Type I error rate for single tests (when the FWER is controlled)
df_TypeI_main_single %>%
pivot_wider(names_from = contrast, values_from = TypeI) %>%
select(Method, N_level, TypeI_cber, everything())
fig_simu_single_b <- df_TypeI_main_single %>%
filter(N_level == 4) %>%
select(-c(N_level)) %>%
mutate(contrast = factor(contrast),
contrast = fct_relevel(contrast, "TypeI_cber"),
contrast = fct_recode(contrast, CBER="TypeI_cber"),
contrast = str_replace_all(contrast, "IV", "level"),
Method = if_else(str_starts(Method, "Main"),
paste0(Method, "\n(uncorrected)"),
paste0(Method, "\n(pairwise only)")),
Method = as_factor(Method),
Method = fct_reorder(Method, TypeI, mean, .desc=TRUE)) %>%
ggplot(aes(x = Method, y = TypeI*100, fill = contrast, group = contrast)) +
geom_col(position = "dodge") +
scale_fill_manual(values = custom_clr) +
geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
labs(x = "One-way ANOVAs with four levels", y = "False Positive Rate (%)", fill = "") +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
theme(legend.position = c(0.5, 0.87),
legend.direction="horizontal",
axis.ticks.length = unit(.15, "cm"),
strip.placement = "outside",
axis.title.x = element_text(margin = margin(t = 5)))
# ggsave(file.path("figures","fig_simu_single_mainpost.png"),
# fig_simu_single_b, width = 12, height = 4.8)
fig_simu_single_b
for (i in unique(df_TypeI_main_single$N_level)) {
nam <- paste("fig_simu_single_b_tmp", i, sep = "")
assign(
nam, df_TypeI_main_single %>%
filter(N_level == i) %>%
select(-c(N_level)) %>%
mutate(contrast = factor(contrast),
contrast = fct_relevel(contrast, "TypeI_cber"),
contrast = fct_recode(contrast, CBER="TypeI_cber"),
contrast = str_replace_all(contrast, "IV", "level"),
Method = if_else(str_starts(Method, "Main"),
paste0(Method, "\n(uncorrected)"),
paste0(Method, "\n(post-hoc only)")),
Method = as_factor(Method),
Method = fct_reorder(Method, TypeI, mean, .desc=TRUE)) %>%
ggplot(aes(x = Method, y = TypeI*100, fill = contrast, group = contrast)) +
geom_col(position = "dodge") +
scale_fill_manual(values = c(custom_clr[c(1,3)], brewer.pal(8,"Set1"),
brewer.pal(6,"Set3"))) +
geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
labs(x = paste0("One-way ANOVAs with ", i, " levels"), y = "False Positive Rate (%)", fill = "") +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
theme(legend.position = c(0.5, 0.88),
legend.direction="horizontal",
axis.ticks.length = unit(.15, "cm"),
strip.placement = "outside")
)
}
fig_simu_single_b_supp <- ggarrange(fig_simu_single_b_tmp3,
fig_simu_single_b_tmp4,
ncol=2) %>%
ggarrange(fig_simu_single_b_tmp5,
fig_simu_single_b_tmp6,
ncol=1)
# ggsave(file.path("figures","fig_simu_single_mainpost_supp.png"),
# fig_simu_single_b_supp, width = 14, height = 14)
fig_simu_single_b_supp